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Abstract. We describe a new data collection scheme for optical diffusion tomography 
in which plane wave illumination is combined with multiple projections in the slab 
imaging geometry. Multiple projection measurements are performed by rotating 
the slab around the sample. The advantage of the proposed method is that the 
measured data can be much more easily fitted into the dynamic range of most 
commonly used detectors. At the same time, multiple projections improve image 
quality by mutually interchanging the depth and transverse directions, and the scanned 
(detection) and integrated (illumination) surfaces. Inversion methods are derived 
for image reconstructions with extremely large data sets. Numerical simulations are 
performed for fixed and rotated slabs. 
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1. Introduction 

Tomographic imaging with diffuse light, often referred to as optical diffusion tomography 
(ODT), is a novel biomedical imaging modality [1,2]. Although ODT was introduced 
more than a decade ago, efforts to bring it into the clinical environment are hampered 
by relatively low quality and spatial resolution of images. Therefore, optimization of 
image reconstruction algorithms for high-resolution ODT is of fundamental importance. 
In this paper we study the image reconstruction problem of ODT by combining three 
novel approaches. First, we employ analytic image reconstruction methods which allows 
the utilization of extremely large data sets [3,4]. Second, we make use of multiple 
projections [5]. Here by multiple projections we mean multiple orientations of the 
measurement apparatus with respect to the medium. Finally, we utilize the recently 
proposed plane wave illumination scheme [6]. Each of these methods provides an 
advantage which is not lost when the techniques are combined. We begin by briefly 
reviewing the approaches to ODT imaging mentioned above. Note that throughout this 
paper we consider the slab imaging geometry which is often used in mammography and 
small-animal imaging [7,8]. In order to obtain multiple projection measurements, a pair 
of parallel plates are rotated around the medium to be imaged which is assumed to be 
stationary and unperturbed. 

There is a direct relationship between the spatial resolution of images and the 
number of data points used for reconstruction [3]. Indeed, the reconstruction of an 
image with N voxels, in principle, requires at least N measurements. In practice, the 
ill-posedness of the image reconstruction problem and the presence of noise require that 
this number be larger than N. Measurements with up to 10 10 data points are feasible 
with CCD camera-based instruments. However, many previous studies of the image 
reconstruction problem in ODT have been limited to relatively small data sets (e.g., 256 
data points in Ref. [9], 900 data points in Ref. [10]). This can be explained by the high 
computational complexity of algebraic image reconstruction algorithms which scales as 
0(N 3 ). To ameliorate this difficulty, we have recently introduced a family of analytic 
image reconstruction algorithms that can utilize extremely large data sets [11-15]. These 
methods allow a dramatic reduction in computational complexity which, in turn, leads 
to a significant improvement of spatial resolution of images. However, these methods 
have certain limitations. 

First, the data collection method described in Ref. [14] requires that measurements 
are taken for source-detector pairs separated by a distance which is much larger than 
the slab thickness. In practice, such measurements are technically difficult to perform. 
Reduction of the required dynamic range of the detectors can be achieved by using plane 
wave illumination [6]. Note that due to the general theoretical reciprocity of sources 
and detectors, plane wave illumination and scanned detection is equivalent to integrated 
detection and scanned narrow beam illumination. However, in a practical situation, the 
different nature of illuminating and detecting devices must be taken into account. For 
the sake of definitiveness, we consider below plane wave illumination and combine it 
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with analytic image reconstruction methods. Note that plane wave illumination requires 
time- or frequency-resolved measurements. However, it can be seen that the number of 
degrees of freedom in the data is still insufficient for unique, simultaneous reconstruction 
of the absorption and diffusion (or reduced scattering) coefficients. This situation is 
similar to the nonuniqueness demonstrated in Ref. [16]. Therefore, we focus here on 
the reconstruction of absorbing inhomogeneities assuming that the diffusion coefficient 
of the medium is constant. Reconstruction of purely absorbing inhomogeneities have 
been employed, for example, in breast imaging [17-20] or blood oxygenation level 
imaging [21,22]. 

Second, it was shown in Ref. [3] that in the slab imaging geometry the depth 
resolution (in the direction perpendicular to the slab) is fundamentally different from 
the transverse resolution (in the direction parallel to the slab surface). The depth 
resolution is much more sensitive to noise and the point-spread functions (PSFs) in the 
depth direction strongly depend on the location of the inhomogeneity. This results in 
image artifacts. In general, the non-uniformity of the PSF can be a serious problem 
if more than one inhomogeneity is present. To correct this situation, we have recently 
proposed multi-projection image reconstruction methods [5,15]. Multiple projections 
render the depth and transverse directions mutually interchangeable. As a result, the 
PSF becomes more uniform and less position-dependent, and also more sharply peaked. 
Note that multiple projections have been used in X-ray imaging for some time. However, 
an important difference between ODT and X-ray computed tomography is that, in the 
first case, tomographic imaging is possible in principle with a single projection while 
in the second case it is not. Perhaps, due to this fact, multiple projections in optical 
tomography have not been investigated until recently, except for the case of ballistic 
propagation without scattering (e.g. [23]), or in conjunction with a modified version 
of X-ray backprojection tomography with phenomenological corrections introduced to 
compensate for scattering [24,25]. In Ref. [15] we have developed a general theoretical 
formalism for inverting measurements obtained from multiple projections. In Ref. [5] 
image reconstruction with two orthogonal projections was numerically implemented. 

In this paper we implement the more general image reconstruction algorithm of 
Ref. [15] for treatment of more than two projections in conjunction with plane wave 
illumination. Note that the plane wave illumination is advantageous when measurement 
are limited by the dynamic range of detectors. If the dynamic range is not an important 
experimental factor, the traditional measurement scheme with point sources and point 
detectors is expected to provide superior image quality. We combine the advantageous 
features of these two approaches with the computational efficiency of the analytic image 
reconstruction methods. 



Multiple Projection Optical Tomography 



4 



2. Theory 

2.1. Single projection 

We assume that propagation of multiply-scattered light in tissue is described by the 
diffusion equation. In addition, we will work in the frequency domain with the sources 
harmonically modulated at the frequency uj and detectors which yield the oscillatory 
part of transmitted intensity. Then the density of electromagnetic energy in the medium 
u(r) obeys the diffusion equation 

- D V 2 u(r) + [a(r) - iu)u(r) = S(r) , (1) 

where a(r) is the position dependent absorption coefficients, S(r) is the source function 
and the D Q is the diffusion coefficient. 

Consider a slab of thickness L with the plane of incidence located at x — —L/2 and 
the detection plane at x — L/2. The medium is located in the region —L/2 < x < L/2. 
If point-like sources and detectors are used (typically, thin optical fibers), the data can 
be expressed as a function (j)(u>, p s , p d ), where p s and pd are two-dimensional vectors 
specifying the location of the sources and detectors, respectively, on the slab surfaces. 
Using the first Born approximation, we linearize the forward model by decomposing 
the absorption function a(r) into a constant background and a small fluctuating part, 
a(r) = ao+5a(r). We seek to reconstruct the values of 8a{r) from the data 4>(u>, p s , pd). 
The usual mathematical formulation of the ODT inverse problem is based on the integral 
equation [26] 

<f>(u, p s , Pd) = J r(w, ps, Pd, r)5a(r)d 3 r , (2) 

where 



tv \ f d 2 q s d 2 q d , 

T(uj,p s ,pd;r) = J K{u,q s ,q d ;x 



(2irY 

x exp [iq s ■ (p - p s ) + iq d ■ (pd ~ p)\ , (3) 

p is the transverse part of the vector r (r = (x, p)) and the form of k(u, q s , qd', x) is 
determined from the boundary conditions on the surfaces of the slab and the expression 
which relates the measurable intensity to the energy density u(r). The derivation of 
(El,© and explicit expressions for k are given in Ref. [3]. Note that the general form 
of P j) .(|3* |) follows from the symmetry of the problem and is independent of the diffusion 
approximation. 

Next, we introduce the plane wave illumination scheme. Instead of using point 
sources located at points p s , we illuminate the slab with a normally incident wide 
homogeneous beam of sufficiently large diameter (compared to transverse dimensions of 
the slab). At the same time we utilize point detectors. This ensures that the new data 
function ip(u,pd) defined by 

ij>(w, p d ) = J <f>(u, p s , Pd)d 2 p s (4) 
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has the same number of degrees of freedom as the unknown 8a(r) (two spatial directions 
and the frequency uj). Thus, the inverse problem is well determined. The integral 
equation (J2J) can now be transformed to 

f d 2 o 

ip(u, p d ) = J j^-^k(uj, 0, q; x) exp[iq ■ (p d - p)}5a(r)d 3 r . (5) 

If ij) is measured for N different modulation frequencies and the sources are placed on 
a square lattice with step size h, Eq. (J3j) can be inverted using the methods described 
in [14]. The SVD pseudo- inverse solution is given by 

f d 2 u ~ 

5a(r) = h 2 j — -^exp(-zn-p)^P*(^n;r)(^|M- 1 (n)|^V( W / ,n) . (6) 

FBZ ^ 71 ) lu^lu' 

Here the vector u is in the first Brillouin zone (FBZ) of the lattice of sources, namely, 
—Tt/h < u y>z < ir/h and 

P(w, u;r) = '^2 K ( w ) 0, u + v; x) exp(iv ■ p) , (7) 

V 

where v are reciprocal lattice vectors of the form v = (2ir / h)(n y y + n z z) . The elements 
of matrix the M(u) are given by 

(u\M{u)\u')=Y J M l {u + v) , (8) 

v 

where 



/■L/2 

(u\Mi(q)\u') = / K(u,0,q;x)K*(u',0,q;x)dx (9) 

J-L/2 



rL/2 

:/s 

(the inverse matrix M~ l (u) must be appropriately regularized [27]) and the Fourier 
transformed data function i/)(u,u) is defined as 

u) = p d ) exp(iu ■ p d ) . (10) 

Pd 

Note that, if 5a is reconstructed only at points which are commensurate with the lattice 
of sources, the factor exp(iv-p) is equal to unity and the function P becomes independent 
of p. Note also that k and Mi can be calculated in terms of elementary functions [15]. 



2.2. Multiple projections 

We now consider inclusion of multiple projections. Let the sources and detectors be 
rotated around the sample as illustrated in Fig. ^ We assume that the rotations do not 
disturb the medium inside the cylinder y/x 2 + y' 2 < L/2 and that the unknown function 
5a vanishes outside the same region. The space inside the slab but outside the above 
cylindrical region is assumed to have the background values of the coefficients cto and D . 
Experimentally, this can be implemented, for example, by rotating an imaging apparatus 
around a sample suspended in matching fluid. We introduce cylindrical coordinates 
r = (R, z, ip) with the z-axis being the axis of rotation. If the data are measured for 
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Figure 1. A sketch of the experimental set up with rotating slab. The axis of 
rotations (Oz) is perpendicular to the plane of the figure and coincides with the axis of 
the cylinder R < L/2 inside which reconstructions are performed. Locations of sources 
and detectors are given in a local reference frame which rotates together with the slab. 

Ng different orientations, where the respective angles 9 n are equally spaced and given 
by n = 2 r w{n — 1)/Ng, n = 1, ... , Ng, the reconstruction formula © can be generalized 
to [15]: 

5a{r) = — — 2^/ — exp[-t(u z z + n<p)\ ^ / du y / duP(uj,u,n;r) 
x(u,u y \M~ l (u z ,n)\uj',u' y )i!(uj f ,u' y ,u z ,n) . (11) 

Here 



P(uj,u,n;r) = J3 q(k>, tt + n + iVflfc; R) exp\i(Ngk(p + ^ z z)] , (12) 

fc=-oo « 
/•27T 

a(u, q,m; R) = / /c(o;, 0, q; iicos ip) exp[i(oyii!sin <p — mip)]dip , (13) 

J 

the elements of the matrix M(u z , n) are given by 



(u,u y \M(u z ,n)\u',Uy) = Y{u,Uy + v y \M l (u z + v z ,n + Ngk)\u}\u' y + v' y ) , 

k=—oov y ,v' y v z 

(14) 

(uj,qy\M 1 (q z ,'m)\u ,q y ) ^ J a(u, q y , q z , to; R)a*(u , q y , q z , to; R)RdR (15) 
and the Fourier-transformed data function is 

u, n) = 53 ^(^, Pd, exp[z(u • p d + n0)] . (16) 
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Note that in (J 16)) we have explicitly included the dependence of the data function on the 
angle of orientation 9. The functions a(u,q,m; R) and (uj,q y \Mi(q z ,m)\uj r ,q y ) can be, 
in general, expressed in terms of modified Bessel functions. The corresponding integrals 
dT3J) and ([To]) are calculated in the Appendix for the case of purely absorbing boundaries. 

A few comments on the reconstruction formula (fTTjl are necessary. First, there is 
an apparent difference between the variables u z ,n and u y ,u>. The first set of variables 
correspond (after Fourier transformation of the data) to the variables z, 9. These are 
the variables with respect to which the unperturbed medium is translationally invariant, 
and they can be referred to as "external" variables. The variables u,u y are "internal" 
variables: they do not correspond to any translational invariance of the system. Second, 
the reconstruction algorithm (JT5j) involves integration over the continuous variables 
u y and u'y and inversion of the operator M(u z ,n) whose matrix elements depend on 
continuous indices. However, if the variables discretized and the corresponding 

integration in (|16|) is replaced by a summation, then M(u z , n) becomes a discrete matrix. 
The resulting reconstruction formula is no longer an SVD pseudo-inverse on the whole 
set of data pa, 9). However, it is a pseudo-inverse solution on the set of the Fourier- 
transformed data ip(uj,Uy,u z ,9) where u y takes only discrete values. Third, it can 
be verified that in the case Ng = 1, the reconstruction formula (fT5j) reduces to ©• 
Fourth, we note that the number of degrees of freedom in the data-function ip is four 
(u,u y ,u z and n). Thus, when the number of rotations is large, it is sufficient to use 
only one or a few values of the variable u y , in which case the inverse problem is still well 
determined. It can be argued that the reconstruction algorithm is then "numerical" 
in one dimension and "analytic" in two. § However, when only a small number of 
projections is taken, we must use a relatively large number of discrete values of u y . By 
doing so, we increase the size of the matrix M whose SVD must be found numerically. 
The inverse solution (JTTj) is then "numerical" in two dimensions and "analytic" in one. 
A similar algorithm (numerical in two dimensions and analytic in one dimension) was 
proposed and implemented in [5], where the image reconstruction area was rectangular 
rather than cylindrical, but only two orthogonal projections were allowed. In contrast, 
the full potential of the image reconstruction algorithm proposed here is realized when 
Ng is large. 

3. Numerical Results 

3. 1 . Single projection 

We have implemented the proposed reconstruction algorithm using computer-generated 
data and the following parameters: the slab thickness was chosen to be the same as the 
cw diffuse wavelength, L = 2tc^J D /a (for most biological tissues, this corresponds to 

§ If the number of rotations and the number of discrete values of u y are both large, it should be possible 
to recover the absorption and scattering coefficients uniquely and simultaneously. This theoretical 
possibility is not discussed in this paper. 
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d=d =0.125L d=d =0.25L d=d =0.375L 




d=d =0.5L d=d =0.625L d=d =0.75L 




Figure 2. Tomographic slices parallel to the slab surface drawn through the medium 
at different depths d (from the plane of scanned detection) with the small absorber 
lying in the center of the field of view at the same depth do = d, and the point- 
spread functions representing depth resolution (a,b). The curves are plotted on the 
same scale (a) and normalized to their own maxima (b). For curves (a,b), the point 
absorber depth is do = 0.251/ (solid line), da = 0.51/ (short dash) and do = 0.75L (long 
dash). 

L ~ 6cm); the lattice step was chosen to be h = L/40 and we have used N = 25 different 
modulation frequencies which range from uo = to uo = 10ao (the maximum frequency 
corresponds to ~ 1.6GHz); the field of view was chosen to be L x L and, finally, we have 
generated forward data for a single point (delta-function) absorber which is located in 
the center of the field of view but at different depths. Absorbing boundary conditions 
were imposed on the surface of the slab. The corresponding expression for the function 
k(u, 0, q; x) is given in the Appendix. 

The results of reconstructions are shown in Fig. |21 The density plots represent 
tomographic slices of the medium drawn at different depths d (the distance from the 
plane of scanned detection) parallel to the slab surfaces. The depth of the absorbing 
inhomogeneity, d , was in each case equal to d; thus the slices represent the depth- 
dependent y — z PSFs. Each density plot has a linear color scale and is normalized 
to its own maximum. As expected, the PSFs become broader when the point absorber 
approaches the illuminated plane. The last two panels (a,b) show the PSFs in the depth 
direction (x) for point absorbers located at d = 0.25L, d = 0.5L and d = 0.75L. 
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Note that the approximate half-widths of these curves are 0.06L, 0.09L and 0.09L, 
respectively. 

The analysis of Fig. El suggests that the PSFs are depth-dependent. Moreover, the 
PSFs have different integral weights. Thus, the point absorbers which are closer to the 
plane of scanned detection result in higher peaks in the reconstructed images. The width 
of the PSFs also depends on depth of the point absorber. This potentially constitutes 
a serious problem for three dimensional tomographic imaging. 

3.2. Multiple projections 

We have implemented numerically the multi-projection image reconstruction formula 
(jlljl . Note that in the multi-projection case there are two choices for graphically 
representing the tomographic slices. In one case, the slices are perpendicular to the 
axis of rotation. The image then is reconstructed in a circle. This choice is convenient 
for studying the radial and angular resolutions. Another possibility is to construct 
cylindrical slices R = -Ri mag e — const, and map them onto rectangles. The image is then 
reconstructed in the rectangular area 27r_R image x (z max — z min ), where z max and z m[n are 
the maximum and minimum values of z, chosen arbitrarily. 

We start with the discussion of circular slices. The results of numerical 
implementation of the reconstruction formula (fTT)|) are shown in Fig. El for four different 
orientations of the slab, namely 9 = 0, tt/2, it, 3ti/2. We have used 23 equally spaced 
values of u y ranging from —ir/h to ir/h and 15 equally spaced modulation frequencies 
ranging from to IOoq] otherwise, the parameters are the same as in Fig. El The 
inhomogeneity was located as specified in the figure legend. The white spots in the 
density plots illustrate the depth PSFs. The graphs (a,b) show the same PSFs in a more 
quantitative way by plotting 5a along the diameter of the cylinder which intersects all 
three inhomogeneities. 

As expected, using four different projections improves the image quality by 
interchanging the source and detector planes, and the depth and transverse directions. 
Moreover, using more projections than four does not change the results substantially, 
as is illustrated in Figs. Eland El However, when a large number of projections is taken, 
the inverse problem becomes well determined even when a relatively small number 
of "internal" degrees of freedom u y is used. This makes the reconstruction formulae 
computationally efficient. Thus, the computation time required for producing data for 
Fig-Elis more than an order of magnitude less than that for Fig. EJ yet the image quality 
is similar. We have verified that three discrete values of u y is also sufficient for Ng = 20 
(taking a single value u y = results in a slight decrease in image quality; data not 
shown) . 

Although the images shown in Figs. E1E1 are similar, the best image quality is, in 
fact, attained in Fig.HJ Here the approximate half-widths of the PSF in the R direction 
are 0.05L for the inhomogeneity located at i?o = 0, 0.04L for the inhomogeneity at 
Rq = 0.25L, ip = 0; and 0.03L for the inhomogeneity at R = 0.375L, ip = ir. These 
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Rn=0 



z = 

R = 0.25L, (p = R = 0.375L, (p =n 




Rn = 



z = . 05L 

R = 0.25L, cp = R = 0.375L, <p =7r 



Ro = 



z=0 . 1L 

R = 0.25L, (p = R = 0.375L, (p =n 





-0.5 0.5 -0.5 0.5 



Figure 3. Circular slices illustrating radial, angular and z resolution. All point 
absorbers are in the z — plane, and the point-spread functions representing depth 
(R) resolution (a,b). The radial and angular coordinates of the point absorber, i? and 
tpo, a re specified in the figure legends. First row of images: slices at z = 0; second row: 
slices at z = 0.05L; third row: slices at z — 0.1L. Images (a-b): reconstruction along 
the diameter that crosses all three inhomogeneities. In (a,b) solid line corresponds to 
i?o = 0.375L and ipo — t, short dash to i?o = and long dash to i?o = 0.25L and 
(^o = Four projections, 15 modulation frequencies and 23 discrete values of u y are 
used. 
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z = 



R = R =0.25L, cp = R =0.375L, q) =n 




z = . 05L 

R = R = 0.25L, cp = R = 0.375L, cp =n 




z = . 1L 

R = R = 0.25L, (p = R = 0.375L, cp =n 




(a) (b) 




-0.5 0.5 -0.5 0.5 



Figure 4. Same as in Fig. [3| but 20 projections are used. 
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z = 



Ro = 




R = 0.25L, <p = R =0.375L, ip =7T 



* M * ' 1 




Ro = 



z = . 05L 

R = 0.25L, cp = R = 0.375L, (p =7T 




Ro = 



z = . 1L 

R = 0.25L, cp = R = 0.375L, cp =n 




(a) 



(b) 




5cf[a.u. 



i \ 



x/i 
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-0.5 0.5 -0.5 0.5 



Figure 5. Same as in Fig. [3] but 40 projections and only three discrete values of • 
are used. 
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values should be compared to the respective values given in the discussion of Fig. 01 
In particular, the inhomogeneity located at h = 0.5L in Fig. El corresponds to the 
inhomogeneity at Rq = in Figs. 0HS1 and is the most "difficult" to reconstruct since it 
is located deep inside the medium. It can be seen that the PSF half-width in the image 
of this particular inhomogeneity is reduced by approximately the factor of 2 due to the 
use of multiple projections. In addition, the relative heights of the maxima of the PSFs 
in Fig. EE do not differ as much as in Fig.EJ This is expected to reduce image artifacts. 

Now we consider the cylindrical slices. From the computational point of view, the 
use of cylindrical slices is a more natural way to display reconstructed images. This 
is evident from the inversion formulae pi |l . p2jl . Indeed, it can be seen that when the 
reconstructed image is rasterized so that the variables z and (p are placed on lattices with 
steps h and 2n/Ng, respectively, the function P(u,u,n; R, z,ip) becomes independent 
of z and (p. Then the dependence of reconstructed images on these two variables is only 
due to the exponent in the integral (jllj) and the reconstruction formula, with respect 
to these two variables, is reduced to a Fourier transform. In Fig. |H1 we have used three 
discrete values of u y with 40 different projections and slices are drawn as described in 
the figure caption. Fig. Efa) illustrates image reconstruction with noiseless data. It can 
be directly compared to slices shown in Fig. El To demonstrate the stability of image 
reconstruction, we have added random Gaussian noise to the data function at the level 
of 1% of the average absolute value of the data. The result is shown in Fig. Efb). As 
is well known, inclusion of noise tends to decrease spatial resolution. It can be seen 
that this effect is stronger for inhomogeneities that are deeper inside the medium. We 
have demonstrated earlier that multi-projection imaging is more stable in the presence 
of noise than the single projection technique [5]. 

4. Summary 

In summary, we have presented a new experimental modality and computationally 
efficient image reconstruction algorithms for optical diffusion tomography employing 
plane wave illumination with multiple projections. Note that due to reciprocity, plane 
wave illumination and scanned detection is equivalent to illumination by a scanned 
narrow beam and integrated detection (e.g., with the use of time-resolved CCD camera). 
The following specific conclusions can be drawn 

• Use of plane wave illumination may be simpler experimentally than the 
traditional approach in which point-like sources and detectors are scanned because 
measurements with a much smaller dynamic range are required. 

• In a single projection experiment, the image quality is relatively good when the 
point absorber is close to the scanned surface and deteriorates as it approaches the 
plane of illumination. This situation should be contrasted with the traditional 
point source/point detector modality [14], where the image quality is low 
for inhomogeneities located in the center of a slab and improves when the 
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R=R =0 . 9 (L/2) R=R =0 . 9 (L/2) 




(a) (b) 

Figure 6. Cylindrical slices illustrating z and ^-resolution for zero noise level (a) and 
for 1% noise-to-signal ratio (b). The point absorbers are located in the z — plane 
at radial depths Rq as indicated. The cylindrical surfaces with radii R = R (directly 
intersecting the inhomogeneity) are shown as projections onto a plane; the length of 
the vertical side of each rectangle is equal to L and of the horizontal side to 2irR. 
Forty projections, 25 modulation frequencies and 9 discrete values of u y are used for 
reconstruction. 
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inhomogeneity approaches either of the imaging surfaces. For a point inhomogeneity 
in the center of a slab, the image quality is slightly better for the traditional (point 
source/point detector) modality (cf. [14]). 

• Rotating the slab around the sample removes many of the deficiencies of the plane 
wave illumination scheme by interchanging the scanned and integrated detection 
surfaces and depth and transverse directions. A minimum of four projections is 
required for such an interchange. 

• When only four rotations are used, a large number of discrete values of the 
"internal" variable u y must be utilized in the reconstruction. Alternatively, a large 
number of projections can be used with a small number of discrete values of u y . 
The second approach is much more computationally efficient but requires more 
complicated measurements. The quality of images is similar in both cases. 

• The plane wave illumination approach allows one to significantly simplify 
reconstruction formulae, both in single- and multiple-projection imaging. 

• If only small number of projections is used (two or four) an alternative approach 
may be used, which is purely numerical in two dimensions and analytic in one 
dimension [5]. For a large number of projections, the algorithm reported here is 
computationally more efficient. 

This work was supported in part by the AFOSR under the grant F41624-02-1-7001 
and by the NIH under grant P41RR0205. 

Appendix: Calculation of the functions a(u,q,m; R) and Mi(q z ,m). 

The function a(u, q, m; R) is defined by (fI3j). To evaluate the integral, we must specify 
the function k(u, 0, q; x). Explicit expressions for k are given in [3] for general boundary 
conditions. In this paper we consider absorbing boundaries for which k is given by the 
expression 



where £* = 3D /c is the transport mean free path, c is the average speed of light in the 



medium, k = y («o — ioj)/Do is the complex diffuse wavenumber, Q = \Jq 2 + k 2 and 
q = (q y ,q z ). Generalization to mixed boundaries of Robin type is straightforward and 
is not discussed here. Then, the expression for a(u,q,m; R) becomes 



x / sxnh[k(L/2 — Rcostp)) sinh[Q(L/2 + Rcostp)] exj>[i(q y Rsxn<p — m(p)]d<p . 




(Al) 





(A2) 
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a(u>, q, m; R) 
x jexp 
— exp 



1 



\D J 4 sinh(fcL) sinh(QL) 
(Q + k)L/2] F m [{Q - k)R, iq y R~\ - exp 

(Q - k)L/2\ F m \(Q + k)R, iq y R] + exp 



-Q + k)L/2 
-Q - k)L/2 



-Q — k)R, iq y R 



-Q + k)R, iq y R 
(A3) 



where 



F m (u,v) 



J exp[u cos <f + v sin if — imipjdip = 2n y 



Vu 2 + V 2 
u + iv 



(A4) 



and I m (x) is the modified Bessel function of the first kind. Note that (jA4|) is well 
defined, including the case v = iu. 

The expressions f)A3|) and (jA4|) define a(q,m;R). Next, we need to calculate the 
matrix elements of M\(q z , m). This integral contains sixteen terms of the form 

4sinh(fcL)sinh(QL)sinh(^)sinh(g-L) ^ ^ + 8 * Q + + "* Q W] 



- Sl k + s 2 Q) 2 - q 2 ] [{-s 3 k> + s 4 Q') 2 - (q' y ) 
(sxk + s 2 Q - q y ){s 3 k' + s 4 Q' - q' y ) 



fL/2 
X / I„ 

Jo 



Rj(- Sl k + s 2 Q) 2 - q 2 I m RJ(-s 3 k' + s 4 Q') 2 - {q' y 



RdR 



(A5) 



where s& = ±1, the sixteen terms correspond to sixteen possible permutations of 
the signs of and the primed variables should be understood as follows: k' = 



a Q — iuj')/D and Q' = \j(q' y ) 2 + q% + (k') 2 . The integral in (jA5|) is evaluated with 
the use of 



xl n (ax)l n (bx)dx 



ij °_ a [al n+1 (ac)l n (bc) - bl n (ac)l n+1 (bc)) , a ^ b 



(A6) 



a = b . 



This completely defines all the functions necessary for implementation of the multi- 
projection reconstruction algorithm. 
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